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ABSTRACT 


Power and fuel temperature excursions in a fission 


realtor 

del*y ed 


in a 


are analyzed using two-temperature model with and without 
neutrons. The governing differential equations are recast 

_ vS tem 

dimensionless form for theoretical studies. The resulting 

cati° n 

of differential equations is analyzed using Hopf bi* arc 
theorem. It is found that depending upon the reactor kinet i cS 
thermal hydraulic parameters, the conditions of Hopf b ifurcat* 
theorem may or may not be satisfied. The parametric equation® 

The 

the bifurcation line or curve are obtained and plotted- 

presence of Hopf bifurcation confirms the existence of P eri 

£ nr- the 

orbits in a finite neighborhood of the operating point, tQl 
appropriate values of the parameters as ascertained from 
bifurcation curve. Stable limit cycles for both with and wit * 1041 
delayed neutrons are calculated numerically, showing eustal** 

nonlinear periodic oscillations in reactor power, fuel and cool &n 

J 4 iod® 

temperatures, and the precursor concentration. The time P ei 

. « I j^ed 

and the amplitudes of the oscillations are also opt -* 1 
numerically. 

It is found that for the conditions of Hopf bl furcati * 2 
theorem to be satisfied, the coolant temperature coefficl® nt 
be negative. Fuel temperature coefficient may be posltl ve ° 
negative. In both cases, stable limit cycles exist. 
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NOMENCLATURE 


reactor power, U 

delayed neutron precursors (in power units), U 
operating or steady state values of P and C, U 
dimensionl elss power [ = (P - P )/P ] 

dimensionless delayed neutron precursors [ = (C-C )/C ] 

o o 

average fuel temperature, °K 

inlet & exit temperatures of the coolant, °K 

average coolant temperature [ = (T e +T^)/2], °K 

operating or steady state values of T r and T , °K 

t c 

dimensionless fuel temperature [= (T f -T°) /(T°-T . ) ] 

dimensionless coolant temperature [= (T -T°)/(T°-T.) ] 

c c c i 

neutron life time, sec. 

fuel and coolant heat capacities, JK 1 
cn £ (T^°-T.)/(X 1), dimensionless 
a (T °-T.)/(X. 1), dimensionless 

C C 1 

ft / (X 1), dimensionless 

C,/C , dimensionless 
f c 

P /[X C (T,°-T °)], dimensionless 
o f f c J 

(T c °-T.)/(T f 0 -T i ) , dimensionless 

fuel and coolant temperature coefficients of 
reactivities, K * 

delayed neutron fraction, dimensionless 
delayed neutron decay constant, sec * 
reactivity, dimensionless 



t dimensionless time ( = X t) 

H f fuel to coolant heat transfer coefficient (UK 

H c Twice the product of the mass flow rate and 

heat of the coolant (UK * ) 


) 

specific 



CHAPTER I 


INTRODUCTION 

In this chapter first of all the subject of reactor dynamics, 
in as much as it is relevant to the present study, will be 
introduced. This will be followed by a review of some papers in 
the related field and finally the present work will be outlined. 

i . 1 PROLOGUE 

To understand the dynamic behavior of a nuclear reactor, or 
more precisely the core of a nuclear reactor, one must know at any 
Instant of time the detailed spatial distribution of the neutrons, 
the rate of power production, the temperature and density of all 
the materials present in the core, the pressure distribution, and 
the flow rates of all the materials present in the fluid or 
gaseous phases. The importance of this knowledge can not be over 
estimated since predictions of the core Integrity and the reactor 
safety as well as control depend on how accurately these 
distributions are known. Reactor dynamic calculations are 
concerned mainly with estimating these distributions with 
reasonable accuracy. It must be pointed out that the amount of 
fissile materials and the amount of fission products in currently 
operating 1000 MU reactors are such that if a severe accident were 
to take place leading to a core meltdown and release of the 
fission products to the environment, grave consequences to the 
population might result. Thus careful consideration is required 



In order to develop methods in reactor dynamics which can predict 
the core behavior with sufficient accuracy and reliability. 

To be able to predict the dynamic behavior of physical 
systems, an appropriate mathematical model has to be developed. 
The model should be a true description of the system or the best 
estimate of the true description. However, the equations 
describing the model have to be of a form and level of complexity 
that can be solved with available resources. These requirements 
contradict each other, and therefore the success of any dynamic 
model in predicting the future of system behavior is strongly 
dependent on the assumptions made in formulating the model and the 
problem to which model is applied. 

The core dynamic process can be represented by the flow chart 
shown. A change in the neutron flux or neutron density 
distribution will cause a change in the reactor power distribution 
which In turn will produce changes in the temperature and pressure 
distributions. As a result, the density of the materials 
Involved, namely, the fuel and the coolant, will change. This 
Includes possible phase transitions and the onset of processes 
such as coolant boiling. As a result, the fluid flow rate, both 
in magnitude and direction will also be affected. Processes such 
as fuel expansion, meltdown, the fuel coolant interaction might 
also be Initiated. Changes in the state variables affect the 
various thermal and hydraulic correlations and the neutron cross 
sections of the materials Involved. These in turn affect the 
neutron distribution and so on. Appropriate justification is 
required to break any of the links in the model. Since 
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A FLOWCHART OF THE RELATIONSHIP BETWEEN PROCESSES 
INVOLVED IN REACTOR DYNAMICS CALCULATIONS • 





the links do exist in the physical reality the development of a 
full and precise model of all the processes presented in Flow 
Chart is not feasible. Such a model contains many complicated 
equations in a rather complex three dimensional geometry. In the 
formulation of a practical solution, therefore, several of the 
processes are usually either neglected or approximated by 
simplified models. For example, in many Instances the effect of 
temperature and density changes on the reactor cross sections can 
be approximated by feedback reactivity coefficients. 

When one Is faced with the problem of solving reactor dynamic 
equations, l.e. to find the core power, temperature (fuel, 
coolant, moderator), flux, enthalpy and flow distributions as a 
function of time, many questions arise and must be answered prior 
to choosing the method of solution and even computer programs to 
be used. Some of these questions are listed below which are 
obviously not complete. 

1. Uhat approximation to the Boltzmann equation should be 
chosen - the neutron transport equation, the diffusion 
equation, or the lumped parameters kinetic equations? 

2. Should the problem be solved in three, two or one spatial 
dimension, or should one use the lumped parameter model 
represented by the point kinetic equations in order to obtain 
the change In time of the power distribution? 

3. Uhat approximation, or rather how many energy groups, should 
one use to represent the neutron energy range? Solution 
range from a single group to few (2 to 6 ) energy groups and 
to large number (30) of energy groups . 



4. Uh&t la the required time span of the solution? Some 
transients last a small fraction of a second; others last 
seconds to minutes before the effects are seen; still others 
may last several hours. Each of these probably requires a 
completely different solution scheme. 

5. How should one account for the significant changes that take 
place in the neutron cross sections used in the diffusion or 
transport approximations? In other words, how does one 
account for the feedback due to temperature and density 
changes in the core? 

Based on the type of reactor under consideration, simplifying 
assumptions and approximations, we end up with a set of equations 
which in general are nonlinear. A linear analysis of such a 
dynamical system is useful only when disturbance is small from the 
operating point since in this case only linear terms are dominant, 
nonlinear terms having insignificant contribution. However, if 
the perturbation is large, that is disturbance is large from the 
operating point, it Is the nonlinear terms which dominate. 
Therefore, In order to have a better understanding of the reactor 
dynamics, we have to look for nonlinear analysis and bifurcation 
theory using state space concepts. Linear analysis shows only 
stable and unstable fixed points, that too, only locally but 
nonlinear analysis may show orbits (such as limit cycles which may 
be stable or unstable) which makes It much more Interesting. 
Apart from limit cycles (periodic orbits) nonlinear analysis may 
also show quasi periodic, aperiodic or chaotic orbits. By 
studying bifurcation we may find out whether such orbits occur In 



a particular dynamical system and if so then for what values of 
parameters or what ranges of parameters. 

From control and safety point of view, if for certain values 
of parameters, a disturbance from steady state point is going to 
take it to limit cycle l.e., ending up with continuous oscillation 
in power, fuel temperature, coolant temperature or precursor 
concentration, then such a phenomenon is certainly worth studying. 
A limit cycle is an Isolated periodic solution of an autonomous 
system represented in the state space by an isolated closed path. 
The neighboring paths are not closed, but spiral into or away 
from the limit cycle. The study of such systems using Hopf 
bifurcation theorem provides us with a clue to existence of limit 
cycles. In Hopf bifurcation, a pair of complex conjugate eigen 
values of the Jacobian matrix crosses the imaginary axis with non 
zero velocity, and makes it possible for the limit cycles to 
exist. 

1.2 TIME DEPENDENT REACTOR KINETICS EQUATIONS 

The neutron kinetics equations, sometimes designated as the 
point kinetics equations, are the most frequently encountered 
equation in reactor dynamics. They can be obtained from basic 
principles as is done in [14, 27, 2d]. They can also be 

precisely derived from the transport or diffusion equation as has 
been done in [26]. 

1.2.1 Space Dependent Dynamic Model 

The time dependent diffusion equation for neutrons in a 
reactor as given in [31] is 



^ = Dv <7 2 N - E v N + S 

St a 


(1-D 


where 


Dvv N dV 


E vN dV 
a 


N(r,t) dV = number of neutrons in a volume element dV at a 
point r at time t; 

= number of neutrons diffusing into dV per unit 
time at time t; 

= number of neutrons absorbed in dV per unit time 
at time t; 

S(r, t) dV = number of neutrons produced in dV per unit 
time at time t; 

= diffusion constant; 

= macroscopic neutron absorption cross section; 

= neutron speed. 

Here two assumptions have been made. They are: 

(a) Neutron curpent density obeys Fick’s law, i.e., 

J = - D v V N (1-2) 

(b) All the coefficients appearing in (1.1) are position 
Independent that is their numerical values are essentially 
averages over the neutron velocity distribution. 
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Writing down diffusion equation and delayed neutron equation 
as given by Hetrick [31], we have 


*N 

St 


sc. 
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= D v v N 


E v N + (1 -ft) k E v N 
a ' <o a 

+ E . X.. C. + S 
i l i o 


(1-3) 
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ft . k E V N - X. C. 
i co a i i 


(1.4) 



where , 

k = Infinite medium reproduction factor; 

OD 

ft ^ - delayed neutron fraction for the i-th precursor; 

ft =£/?,; 

i 1 

X. = decay constant of i-th precursor; 

C^(r,t) = Volumetric density of the i-th of precursor; 

S^(r,t) = term representing sources of neutrons different 

from the fission processes. 

Equations (1.3) and (1.4) constitute a set of coupled partial 
differential equations. They together form an infinite 
dimensional dynamical system. Fortunately under certain 

assumptions this system of infinite dimensions can be reduced to a 
finite dimensional system, i.e., a set of ordinary differential 
equations. This can be achieved if we use the lumped parameter or 
point reactor model. 

1.2.2 The Point Reactor Model 

The fundamental assumption underlying this model la that all 
except the fundamental spatial harmonics of the flux, following a 
flux disturbance, have been attenuated in a very short time (a few 
micro seconds) compared to times of reactor kinetics of Interest 
as pointed out by Ash [14]. In other words the density of 
neutrons in the reactor is assumed to rise or fall as a whole 
following a change in reactivity. This means all spatial 
dependences factor away, leaving only time as the Independent 


variable . 



Substituting as in [30], 


N(r,t) = £(r) n(t); C.(r,t) = g.(r) (^(t) (1.5) 

in Eqs. (1.3) - (1.4) and forcing the condition that space 

dependence from the resulting equations must be removed it is 

2 

required that f/g^, V £/£ and S q / f be independent of position. 
Assuming 

g. = f; V 2 £ + B 2 £ = 0 

and letting 

q(t ) = S o (r, t ) / f (r ) 

we end up with the point reactor model as 


where, 


P 

1 

q 


dn 

dt 

dC. 

l 

dt 


€-H n + S X. C. + q 

l i 1 1 


■-i n - X. C. 

1 li 


fundamental buckling, 
reactivity, 

neutron generation time, 
effective source strength. 


( 1 . 6 ) 

(1.7) 


Eqs. (1.6) - (1.7) can be further simplified by assuming 

(a) Reactor is operating at high power, 

(b) All the delayed neutrons have a single mean life time. 


First condition may be used to drop q in (1.6) [30] and the 
second condition implies ft ^ - ft. Then we get the simplified 
equations . 



( 1 . 8 ) 


dn 

dt 

dC 

dt 


_ P- ft 


n + X C 


_ ft 


n - X C 


(1.9) 


Eqs. (1.8)-(1.9) a re the simplified point kinetic model 
which is a system of two coupled differential equations. All the 
assumptions made in arriving at these equations have been clearly 
stated at appropriate places. 


1.3 LITERATURE REVIEW 

In one of the Important papers, Popov [ 2 S ] has shown the 
effect of delayed neutrons on the stability of the reactor. He 
concludes that if stability can be proved by a particular type of 
Liapunov function without considering the effect of delayed 
neutrons, then the system is asymptotically stable when the effect 
of the delayed neutrons Is taken into account. Baran and Meyer 
[22] point out that Liapunov function considered by Popov was of a 
very particular type, and therefore, one cannot assume in general 
that delayed neutrons have a stabilizing effect. Further, he 
gives an example to show that stability without delayed neutrons 
does not imply stability with delayed neutrons. He has also 
derived a new criterion for stability in the large. Gyftopoulos 
et al [32] have shown by considering Lyapunov’s second method that 
If a nuclear reactor is stable for arbitrary small perturbation 
without delayed neutrons, it is also stable for arbitrary small 
perturbation with delayed neutrons. Smets [23] studies linearized 
reactor dynamics equations and shows that delayed neutrons do not 
always improve the stability of nuclear reactors and that a 
reactor may be unstable although it 1s stable when delayed 



neutrons are neglected. Devooght et al [7] discuss four methods 
of analysis of nonlinear stability that can be used when the 
reactors are unstable above a certain equilibrium power. He 
further asserts that topological methods and Liapunov's second 
method may not be applicable for complex systems. For complex 
systems he describes two methods one by Alzermann and second due 
to Rosen. Akcasu and Shotkln [2] studied bounded periodic 
behavior of the reactor power for the cases when the equilibrium 
power is greater than the critical power level in the absence of 
delayed neutrons. For such reactors he also derives expressions 
for amplitude and frequency of the limit cycles. The technique 
they have employed is essentially perturbation method. They have 
also been able to derive a describing function which plays 
Important role in controls. Shotkin et al [3] have considered 
reactors with feedback linearly proportional to the power level 
and shown that consideration of delayed neutrons permit the 
existence of unstable limit cycles. Technique adopted by them is 
essentially the same as [2]. 

Vreeke and Sandqulst [8] present a generalized approach to 
phase space analysis of reactor kinetics equations. Their method 
of approach employs a digital computer and a graphics terminal to 
construct rapidly and accurately two-dimensional and perspective 
three-dimensional phase space portraits over the region of 
interest. In this paper they have considered phase space 
solutions of second and third order systems of equations. The 
method is demonstrated on three different nonlinear problem of 
Interest to nuclear reactor kinetics. A second order problem with 
temperature dependent reactivity and two third-order problems with 



reactivity a function of two effective temperatures and Xe 

concentration are considered. Schmidt and Hetrick [4] use Fourier 
Series methods for determining the frequency and fundamental 
amplitudes of the oscillations in reactivity and power level. 
They also examine the system for the effect of delayed neutrons on 
oscillatory solutions and conclude that when delayed neutrons are 
included, unstable limit cycles as well as stable limit cycles are 
found. Poddar et al [5] have examined oscillatory behavior in 
nuclear reactors ith two temperature feedback using perturbation 
method, but without delayed neutrons. They establish necessary 
and sufficient conditions for asymptotic stability under 

oscillatory conditions. Poddar and Trasi [9] have pointed out 
some discrepancy in [34] and have modified conclusions of [34] . 

By using both Uelton’s criterion as well as Liapunov’s second 

method, the system considered in [34] can only be shown to be 

asymptotically stable in the small and in the large when the 
delayed temperature coefficient lies between certain limits, which 
Include a range of negative and positive values. Gupta and Trasi 
[10] have, for a Newtonian feedback and a step input of 
reactivity, solved point reactor kinetics equations using 

perturbation techniques in the presence of delayed neutrons. They 
have shown that asymptotic power level is Independent of 

perturbation modes higher than the first and that the change in 
the steady power level is inversely proportional to the magnitude 
of the reactivity coefficient of fuel temperature. Uard and Lee 
[6] use singular perturbation method for the analysis of large 
power oscillations in nuclear reactors. Such large oscillations 
occur as the system reactivity becomes significantly positive and ! 



then negative. This paper introduces a potential method for the 
direct limit cycle solution, in phase space. Keener and Cohen [1] 
have shown existence of limit cycle in two-temperature feedback 
model without delayed neutrons using perturbation techniques. 
Ashok Kumar [17] has analysed the two-temperatur e feedback model 
with delayed neutrons. He uses a perturbation approach for 
linearisation and establishes the conditions for oscillatory and 
non oscillatory transient behavior as well as sustained 
oscillations, and verified the results of this linear analysis by 
numerical computation. 

In most of the works just described presence of delayed 
neutrons have been considered only for stability analysis. Limit 
cycle analysis has been done by different techniques neglecting 
the presence of delayed neutrons. Manmohan Pandey [24] has shown 
the existence of stable limit cycles in two temperature feedback 
model neglecting the presence of delayed neutrons by using Hopf 
bifurcation theorem. However, he considered only positive fuel 
temperature coefficient of reactivity. 

1.4 PRESENT WORK 

The purpose of the present study Is to analyze the nonlinear 
system of equations governing reactor power and temperatures, 
without resorting to linearization or perturbation methods. This 
ensures that the predicted behavior will be observed In the actual 
nonlinear system and will be valid for arbitrary disturbances. 

The model used is the two-temperature feedback model, with 
and without delayed neutrons. For each of the models, first. 


conditions for sustained oscillations will be established. 


Then 



It will be shown that all the conditions o£ Hopf bifurcation are 
satisfied and hence either a stable limit cycle or an unstable 
limit cycle or a nonlinear center exists in the neighborhood of 
bifurcation line/curve. Limit cycles will be located by numerical 
experiments for different values of feedback coefficients. 

In the last chapter, results obtained will be interpreted and 
further investigations required will be pointed out. 



CHAPTER II 


FORMULATION OF DYNAMICAL SYSTEM FOR FISSION REACTORS 


In this chapter, the point reactor model with two energy 
balances and reactivity feedbacks is formulated. These equations 
are put in a convenient form by suitable translation and 
normalization of the dependent and Independent variables. All the 
variables are then non-dimensional. 


£.1 THE LUMPED PARAMETER DYNAMICAL MODEL 

In section 1.2.2, we obtained Eqs . (1.8) and (1.9) for the 
point reactor model with one group of delayed neutrons. Urlting 
reactor power, P instead of neutron flux density n, and using 
power units for C, we can rewrite the equations as 


dP 

dt 


P + XC 


( 2 . 1 ) 


dC 

dt 



( 2 . 2 ) 


Ue shall now Incorporate equations obtained from energy 
balance to Eqs. (2.1)~(2.2) and include their feedback effects in 
the first equation. 


2. 1 . 1 Energy Balances 

The lumped parameter thermal model is based on two energy 
balances: one for the energy contained in all’ the fuel elements of 



the reactor core, and another for the energy stored In the coolant 
within the core volume. 


The energy balance for the fuel elements is given by 


C 


f 


dT f 

dt 


P - H f (T f - T c ) 


where, 

C i 


T , , T 
f ’ c 


heat capacity of the fuel elements in the reactor 
core (J K * ) , 

fuel to coolant heat transfer coefficient (U K * ) , 
average temperatures of the fuel and the coolant 
respectively. 


Dividing all the terms by , we can write the equation as 
dT < 

= KP - r (T f - t c ) 


dt 

where. 


(2.3) 


K = 1/C f and r = H f / C f 


The energy balance for the coolant within the core is given 
by 


C c inr - H f < T f - T C J 

and H = 2 m C 

c c c 

where. 


H 


(T - T.) 
v c i' 


C c = heat capacity of the coolant within the reactor core 

(J K' 1 ), 

m = coolant mass flow rate, 
c ’ 

T^ = coolant inlet temperature, 

= twice the product of the mass flow rate and specific 
heat of the coolant (U K *). 



t 

Dividing all the terms by C , we can rewrite the equation as 
dT 

3 ^ = >-= CT f - T e ) - - (T c - V (2.43 

where c = C,/C , s = H /C and y has been defined in Eq. 

f c c c 

(2.3). 

2.1.2 Reactivity Feedback Effects 

The reactivity p( t) that derives a transient is the net 
effect of contributions arising from the external disturbances to 
which the transient is due and the thermal hydraulic feedback 
effects. Therefore, it can be written as 


p( t) = p. (t) + p fb (t) (2.5) 

where p^(t) is the reactivity from external disturbances and 
p£ b (t) arises due to thermal hydraulic feedback effects. Assuming 
linear feedback, we will take 

Ofb (t) ■ “f (I f - T f“> + “c (T c - T c R > 
where, 

ot c = fuel and coolant temperature coefficients of 
reactivity, respectively, and 

R R 

T^ , T c = some reference values of T^, T c respectively. 


For convenience we will choose these reference temperatures as the 
steady state operating temperatures T^° and T c ° . Then p ^ turns 
out to be zero for Impulse or step reactivity Inputs. Therefore, 
for linear feedback model Eq. (2.5) takes the form 



P 


( 2 . 6 ) 



Here p is time independent. In actual reactors p is a function of 
time. The form of the function may be very complex and even 
unknown . 

For higher order feedback, we can write 

P = « f (T f -T f °) + « c (T c -T c °) + F (T f -T f °, T^ 0 ) (2.7) 

where F (T..-T, 0 , T -T °) is a nonlinear function with no 

I t c c 

constant or linear terms. 

2.2 REDUCTION TO NON-DIMENSIONAL FORM 

Equations (2.1) - (2.4) together with Eq. (2.6) represent the 
lumped parameter dynamic model with linear reactivity feedback. It 
has two fixed points; the shutdown point (0, 0, T^ , T^ ) and the 
operating point (P C , T °, T °) , where P and C are steady 

« w Jl Lt w U 

state operating values of P and C, respectively, of the system 
after disturbance. Now, we shall non-dimensionallze these 
equations with appropriate normalization and translation. 

2.2.1 Normalization of State Variables 

In Choosing normalization and translation of state variables 
following points were kept in mind: 

(1) Resulting transformed state variables should be 
dimensionless , 

(2) The normalization should not be done with respect to 
parameters passing through zero, 

(3) The number of parameters in the final equations should be 


minimum 



( 4 ) 


The operating point should be brought to origin, 

(5) the shutdown point should be independent of the parameters, 
if possible, 

(6) The transformed variables and the parameters should allow 
some direct physical interpretation, as for as possible. 


The normalization and translation used are 


IP 



<c = 


¥ 


f 


¥ = 
c 





( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


2. 2. 2 Time Normalization 

In deciding time normalization, the crlterlons were: 

(1) A dimensionless time should be used, 

(2) Number of parameters should be minimum in the equations, 

(3) The orders of magnitude of the final parameters should be 
reasonable , 

(4) The time scale should not become too Inflated, so as to 
require too large integration steps. 

Time normalization used was: 


T 


xt 


( 2 . 12 ) 



2. 2.3 Thd& Dimensionlesfi Equations 

Incorporation o£ these transformations in Eqs . (2-1) to (2-4) 

and Eqs. (2.6) we obtain the following system of equations: 


dIP 

dr 

dC 

dT 


d¥ 


f 


dr 


d¥ 

_ c 

dr 


- HP + bC + a f TT f + a c V c + a f «P¥ f + ^IPT c 
IP - c 

p(l - S) IP - p ¥ f + p & TT c 


(2.13) 

(2.14) 

(2.15) 

(2.16) 


A total o£ six parameters have appeared in the system of 
equations. They are 


b = ft /\1 , a f = « f (T f ° - T.)/X1, 


«c < T 

° - T.)/A1, 
C 1 

p = H f /C f A, 




(T c ° 

" T ) 





(T f ° 

1 . c 

- V 

= V C c = 

1/K C c 



these 

parameters 

as well as 

state 

variables 

are 


dimens ionl ess . 


In terms of dimensionless IP, C, ¥ ^ and ¥ c> the reactor 
operating or steady state point is (0, 0, 0, 0) and the shutdown 
point is (-1, -1, -1, -1). This suggests that the values of state 
variables below -1 would be unrealistic. 

2.3 TWO-TEMPERATURE FEEDBACK MODEL WITH DELAYED NEUTRONS 

The model used in the previous section was based on 
two-temperature feedback via two energy balances. Delayed 



neutrons were Included and feedback assumed was linear. 

This model already given in previous section will be 
reproduced once again. 

-bCP + bC + a, ¥, + a IT + a. IF¥ . + a 1P¥ (2.17) 

ff C C ff c c 

IP - C (2.18) 

p(l - 6) P - pl f + p6I c (2.19) 

(¥ f - ¥ c ) (2.20) 

This system of four coupled nonlinear differential equations 
continuous dynamical system with six parameters and having at 
least two fixed points (0, 0, 0, 0) and (-1, -1, -1, -1)- The 
fixed points are found by putting derivatives in Eqs . (2.1) 
(2.4) equal to zero. 


dtP 

dr 


dC 

dr 



dT 


d¥ 

c 

dr 


2.4 TWO-TEMPERATURE FEEDBACK MODEL WITHOUT DELAYED NEUTRONS 

If we neglect the delayed neutrons, the parameter b is zero, 
and Eq. (2.18) is no longer required. Thus the governing system 
of Eqs. (2.17) - (2.20) become 


dJP 

- 3 — = a, ¥ , + a ¥ + a, QP¥ , + a IP¥ 

dr ff cc ff cc 


d¥ f 

dr 


p(l * B) IP - p ¥ f + p & ¥ c 


d¥ 

< 

dr 


(¥ f - 


( 2 . 21 ) 

( 2 . 22 ) 


pc 

& 


(2.23) 



This set of differential equations forms a three-dimensional 
continuous dynamical system with five parameters having the two 
obvious fixed points (0, 0, 0) and (-1, -1, -1). 



CHAPTER III 


INVESTIGATION OF THE DYNAMICAL SYSTEMS 


In the last chapter dynamical systems o£ three and four 
dimensions were obtained. Now the behavior of these systems will 
be investigated by analytic as well as numerical methods. The 
main objective will be to locate limit cycles if any. The main 
concept of analysis, which is essentially Hopf bifurcation 
theorem, has been taken from Refs. [15, 20, IS]. 


3.1 TWO TEMPERATURE FEEDBACK MODEL WITHOUT DELAYED NEUTRONS 
This Is a three dimensional system of equations. 


cfiP 

= a IT + a ¥ + a 2F¥ + a EP¥ 

ax £ £ cc f £ cc 


(3.1) 


d¥ 

^ = p(i - 0) IP - P ¥ f + p e ¥ c 


(3.2) 


d¥ 

dT 


- ■ f 2 CW, - 


(3.3) 


It has two fixed points (0, 0, 0) and (-1, -1, -1). These 

are the operating and the shutdown points, respectively. 

Local Behavior at the Operating Point CO, O, O) 

The Jacobian matrix at the operating point is 

0 a, a 

f c 

p(l-©) -p p& 

0 pc/d -pc/O 



Characteristic equation of this matrix is 


X 3 + p(l +c/0) X 2 + (pc/O - a f ) p (1-S) X 

- (a + a ) p( 1 - & ) pc/B = 0 (3.4) 

r c 

Ue write (3.4) as 

X 3 + eX 2 + fX + g = 0 (3.5) 

where, 

e = p (l + c/e), f = (pc/e - a f ) p ( l -e ) , 
g = - (a f + a c ) p(l - e) pc/e (3.6) 

Since sustained oscillation is a precondition for the 
existence of limit cycles, we shall examine the conditions of 
asymptotic stability or sustained oscillations. This can be done 
with the help of Routh - Hurwitz criterion [31]. 

3.1.1 Routh-Hurwitz Criterion 

The characteristic equation is 
3 2 

X +eX +fX + g = 0 

f 
£ 

(3.7) 

From Routh's table it is obvious that all the roots of Eq. 
(3.5) will have negative real parts iff the following conditions 
are fulfilled simultaneously: 


Routh* s Tabul all on 
% 3 


X 

X 


e 

ef-g/e 

£ 



e > 0, 


ef-g/e > 0, 


g > 0 


Sustained oscillations occur when there exists atleast one 
pair of conjugate imaginary roots; the remaining root has to be 
negative real. Ue get information about sustained oscillations 
by looking at the Routh table given above and by the following 
analysis . 


Uhen ef-g = 0, the X - row in (3.7) consists only of zeros, 

so the entry in X° - row will become indeterminate. Uhen a row of 

zeros occurs, there exists an even polynomial divisor (also called 

'auxiliary polynomial’) of the original polynomial. The 

coefficients of this divisor polynomial are given by the previous 

2 

non-zero row of the array. In our case the coefficients in X 
row are e and g, so the divisor is: 

eX 2 + g (3.8) 

This divisor represents a conjugate pair of imaginary roots iff e 
and g have the same sign. 

L'rt’ 

Uhen Routh- Hurwitz array terminates per maturely. Ue 

A 

complete the table by replacing the row of zeros with the 
coefficients of the derivative with respect to X of the divisor: 


dX 


[®x + g] 


The resulting Routh tabulation takes the form 
„ 3 


X 

X 


1 
e 
2 e 


f 

£ 


(3.9) 


£ 


(3.10) 



From the above analysis and the Routh tabulation given in 
(3.10), it follows that the sustained oscillations occur when the 
following conditions are fulfilled simultaneously. 

e > 0, ef - g = 0, g > 0 (3.11) 

3.1.2 Routh-Hurwitz Criterion Applied to the Present System 

From Eq . (3.6), we see that e is always greater than zero, 

since for any reactor p, c, & all are positive definite. Also & 
is always less than unity which means g > 0 is equivalent to 

- ( a £ + a c ) > 0 

or (a £ + a c ) < 0 (3.12) 

Substituting expressions for e, f and g in ef-g = 0 and 

simplification leads to 

?*£-<>„ = P <1 + f> (3-13) 

Eqs . (3.12) and (3.13) are plotted in Fig. 2 for parameter 
values p = 2, C = 1/2, & = 1/20 which represents a typical thermal 
reactor . 

3.1.3 Hopf Bifurcation at Operating Point 

Eqs. (3.12) and (3.13) are necessary and sufficient 

conditions for the roots of Eq . (3.5) having a pair of pure 

imaginary roots. Now treating p, c, fi, a £ to be fixed it is not 

difficult to show that the complex conjugate eigenvalues cross 

the imaginary axis from left to right and the real eigenvalue 

* 

remains in the left-half plane as a £ increases beyond a £ , where 
* 

&£ represents a point on the bifurcation line satisfying Eqs. 
(3.12) and (3.13). To show this we start by assuming that two of 



the eigenvalues are a complex conjugate pair, R ± il, and the 
third o£ course is necessarily real, say, S. Then, by writing the 
characteristic equation, it is straight forward to see that the 
following must hold [24]: 

2R + S = - p (1 + c/0) 

R 2 + I 2 + 2RS = p(l-0) (pc/S - a f ) 

(R 2 + I 2 ) S = p 2 (a. + a ) (1 - 0) cj& (3.14) 

t c 

In order to see the movement of the eigenvalues across the 
Imaginary axis, we differentiate Eqs. (3.14) with respect to a 
parameter, say a^, considering the components R, I and S of the 
roots as variables, and keeping all other parameters fixed. The 
resulting equations can be solved for the derivative of R, and 
considering the fact that c > 0, p > 0, and 1 > & > 0, it is not 
difficult to show [24] that 

> 0 (3.15) 

R=0 

Therefore from Hopf bifurcation theorem it follows that there 

exists either a stable limit cycle for a„ > a * or an unstable 

ft 

limit cycle tor a. < a * or a nonlinear center at a^ ■ a * and all 

ft ft 

the three cases are mutually exclusive. 

Numerical Experiments 

Numerical experiments were performed for p = 2 , c = 0.5, & = 

0.05 for three pairs of values of a^ and a .^ . Pandey [24] reports 
stable limit cycles for negative a £ and positive a^ . Ue 
concentrated on the case where both a c and a^ are negative since 




negative are 


from reactor safety point of view both a c and a^ 
important . 

Limit cycles have been plotted in figures, as Indicated in 
the Table 1. A point within the cycle spirals out tending towards 
limit cycle. A point from outside spirals Inwards towards the 
limit cycle. 

Results of numerical experiments are summarized in the 
Table 1. 


3.2 TWO TEMPERATURE FEEDBACK MODEL WITH DELAYED NEUTRONS 
This is a four-dimensional system of equations 

^ = -WP + bC + a, ¥, + a ¥ + a. IP¥ , + a JP¥ 

dr ffccffcc 


dC 

dr 

d¥, 

dr 

d¥ 

< 

<±T 


= IP - C 


p(l - 6) IP - p T f + p d ¥ c 


P£ (T - ¥ ) 

& K f c J 


(3.16) 


having two fixed points (0, 0, 0, 0) and (-1, -1, -1, -1). 


Local Behavior of the Operating Point CO, O, O, 03 
The Jacobian matrix at the operating point is 


-b 

1 

P(l-S) 

0 


b 

-1 

0 

0 


0 

-p 

PC/S 


0 

p& 

-PC IB 


(3.17) 



Characteristic equation of this matrix is 


X 4 + eX ^ + f X ^ + gX + h ~ 0 (3.18) 

where, 

e = c|£+p+b+l) 

f = p(l -■$) (g^ ~ * £ ) + (b + D P Cl + f) 

g = P£ (b+1 ) p(l-e) -a £ (|H + i) p (l-©) - a c p(l-d) 

h = - |£ p(l-e) Ca f + a c ) (3.19) 

For sustained oscillations like three-dimensional case, Eq . (3.18) 

should also have at least a pair of pure imaginary roots. The 
other pair can either be complex conjugates with negative real 
parts or negative reals. The necessary and sufficient conditions 
for a pair of roots to be complex conjugate can be determined 
using Routh-Hurwitz criterion. 


3. 2 . 1 Rourth-Hurwitz Criterion 

The characteristic equation is 

M O <5 

X + eX + fX + gX + h = 0 

Routh* s Tabulation 

1 f h 

e 8 

ef-g/e h 

g-e 2 h/(ef-g) 
h 

From Routh's table it is obvious that all the roots of Eq. (3.18) 
will have negative real parts iff the following conditions are 
fulfilled simultaneously. 




e > 0; ef-g > g; g - > 0 and h > 0 (3.20) 

Ue can get information about sustained oscillations by 

looking at the Routh table given above and the following analysis. 

2 

When g(ef-g) -e h = 0 , the X - row consists solely of zeros, 

o 

so the entry in the X - row would be indeterminate. When a row 

of zeros occurs, there exists an even polynomial divisor (also 

called 'auxiliary polynomial') of the original polynomial. The 

coefficients of this divisor polynomial are given by the previous 

nonzero row of the array. In the present problem, the 

2 

coefficients in the X -row are ef-g/e and h, so the divisor is: 

ef ~ ? X 2 + h (3.21) 

e 

This divisor represents a conjugate pair of imaginary roots 
iff (ef-g) and h have the same sign. 

To compute the Routh-Hurwi tz array after a p, armature 
termination, we replace the row of zeros with the coefficients of 
the derivatives with respect to X of the divisor: 

d , ef— g . 2 „ 

dX C i - X + h 3 (3.22) 

The resultant Routh tabulation is as follows: 

1 f h 

e g 

ef-g/e h 

2(ef-g)/e 



h 


(3 . 23) 



From the above analysis and the Routh tabulation given 
(3.23), it follows that the sustained oscillations occur when the 
f ol lowing conditions are fulfilled simultaneously: 

g > 0, h > 0 and g(ef-g)-e 2 h = 0 (3.24) 

Using the reactor parameters p = 2.0, c = 0.5, & - 0.05, b = 
88.0, from Eqs . (3.19) and (3.20) we get 

3382 - 39.9 a, - 38 a >0 (3.25) 

f c 

- 38 (a f + a c ) > 0 (3.26) 

a 2 - 0.21 a 2 + 0.73 a, a - 1620 a^ 
f c f c f 

- 1440.5 a + 135985.8 = 0 (3.27) 

c ^ ' 

Eq. (3.27) represents a two-sheet hyperbola which has been 
plotted in Fig. 8. Imposition of conditions (3.25) and (3.26) on 
it gives the bifurcation curve. Three possible zones are shown in 
the figure. 

For all possible values of a^ and a^ on bifurcation curve, 

characteristic Eq. (3.18) has a pair of pure imaginary roots. The 

other two turn out to be pure real negatives. Now in order that 

all the conditions of Hopf bifurcation theorem Is satisfied what 

remains to be shown Is that complex conjugate eigenvalues cross 

the imaginary axis from left to right as one of the parameters (a £ 

or a c ) is increased. To show this analytically involves 

unnecessary algebra. Ue can show the crossing numerically. For 

* 

this we evaluate the roots of characteristic equation at a^ , a f > 

* * 

a^ and a^ < a^ for a number of points and observe that indeed 



crossing is occurring. Crossing of eigenvalues as well as Movement 
of negative real roots are shown in Fig. 7. 

Numerical results for a few points are given in the 
Appendix C. 

Therefore, from Hopf bifurcation theorem it follows that 
there exists a limit cycle in a sufficiently small neighbourhood 
of bifurcation curve. Limit cycle is stable for a^ > 
unstable for a^ < a^.* and for a ■ a^.* there is a nonlinear center 
and all these three cases are mutually exclusive. 

Numerical Experiment 

Numerical experiments were done for p = 2, c = 0.5, & = 0.05, 
b = SS.O for five pairs of values of a^ and a c . Limit cycle 
obtained was stable in all the five cases. 

Results of numerical experiments are summarized in Table 2. 

Starting from any point in four dimension leads us 

asymptotically to the stable limit cycles provided we are not 

too far from the operating point. If we start from a point within 

the limit cycle as seen on ¥ - IP projection it tends to the limit 

cycle. Similarly, if we start from outside as seen in - IP 

projection, the path may appear to cut the limit cycle but in 

reality, that is in four dimension, this is not the case. 

dX 

Otherwise, our system -r— = f (X) will take more than two values 

UT — — 

for the same values of the coordinates. Thus violating the 
single valuedness of the derivative. 



The results obtained in case of four-dimensional case 


shows 


some remarkable differences from the three-dimensional case. 
Firstly in the case where delayed neutrons have been taken into 
consideration has smaller time period than when delayed neutrons 
are neglected. Also amplitude becomes too large in case of model 
with delayed neutrons. 

3.3 NUMERICAL METHOD EMPLOYED FOR SOLVING THE SYSTEM OF 
DIFFERENTIAL EQUATIONS 

The Jacobian matrix for our dynamical system have the 
eigenvalues which are very widely separated. This implies, as 
pointed out in introduction that some transients may last for a 
few seconds and some may last for several hours. Such systems 
cannot be solved quite accurately by methods such as Runge-Kutta 
or Adam-Bashf orth which is based on multiple pr edictor-corr ector 
method . 

To solve our system, the subroutine DG2EBF from NAG library 
was employed on DEC-1090. Subroutine D02EBF is based on Gear 
method. In case of three-dimensional system tolerance used was 
0.1E - 05 and in case of four-dimensional system It was 0.1E-07 
which gives quite accurate solution as was seen by using tolerance 
smaller than these. 

The search for limit cycle was initiated by using a point 
close to operating point as Initial condition. The last step of 
the iteration was again used as the initial condition and it was 
repeated till limit cycle was finally located. Starting from any 
other point asymptotically converges to the limit cycle. To have 
better accuracy the step size were taken small. In case of 



three-dimensional system step size was 0.01 and in 
four-dimensional case, step size was 0.001. 

NAG subroutine C02AEF was used to show numerically the 
crossing of eigenvalues in four dimensional system. The 

subroutine C02AEF calculates all the roots of a real polynomial of 
order upto hundred using the method of Grants and Hitchins. 

The next chapter briefly summarizes the findings of this 


work. 



CHAPTER IV 


CONCLUSIONS 

4..1 SUMMARY OF INVESTIGATION 

Ue started with two temperature model, l.e. two energy 
balances together with reactor kinetics equations to get a set of 
differential equations. Further these differential equations were 
coupled taking feedback proportional to fuel and coolant 

temperatures. This feedback essentially arises because of 

thermal-hydraulic effects. The resulting set of equations were 
non-dlmensionallzed (since it is convenient to analyze 
non-dimensional equations) keeping into mind that number of 
parameters is as small as possible, new state variables have some 
direct physical Interpretations. 

After this, models with and without delayed neutrons were 

analyzed for Hopf bifurcation. Using Routh-Hurwitz criterion, the 

necessary and sufficient conditions for the characteristic 

equations to have a pair of pure imaginary roots were obtained. 

Parameters a. and a were treated as variables where as other 
. f c 

parameters as fixed. Further, it was shown that with respect to 
parameter a f crossing of imaginary eigenvalues is occurring while 
others remain on right hand plane. For model without delayed 
neutrons crossing was shown analytically . For model with delayed 
neutrons crossing was shown numerically. Then using Hopf 





bifurcation theorem the existence of limit cycles were predicted. 
Numerical experiments were done to locate these cycles and were 
found to be stable in both the cases that is with and without 
delayed neutrons. 

4.3 FURTHER INVESTIGATIONS REQUIRED 

Uhile studying the dynamics of a thermal reactor we have 
concentrated on the core of the reactor. Ue have neglected the 
interactions of controls and plant with the core. Further studies 
in this direction can be done taking these Interactions into 
account . 

4.3 EPILOGUE 

In the present work two models, one without delayed neutrons 
and the other with delayed neutrons, were studied using the 
concept of Hopf bifurcation. For the conditions of the theorem to 
be satisfied, and therefore the limit cycles to exist, it was 
found that coolant feedback coefficient has to be negative. Fuel 
feedback coefficient may be negative as well as positive. For 
both the models stable limit cycles were computed numerically. 



o / 


APPENDIX A 


1.1 STIFF SYSTEMS 

"Stiffness" is a property of a mathemtical problem (not of 
the numerical solution method). Let us consider the system of 
differential equations which are first order, nonlinear and 
ordinary. 

d v 

as - f(x - y) ■ 0 y <a) ' y o’ x x [a - *>>• 

The above initial value problem is said to be stiff in an interval 

I <=. [a, bj if, for x £ I, 

(1) Re(X.) <0 (i = 1, 2, ... s); and 

(2) S(x) = max Re ( -X . ) / min Re ( -X . ) >> 0 

i = l,s 1 1=1, s 1 

where the X^ are the eigen values of df/dy evaluated on the 

solution y(x) at x. 

The ratio S(x) may be termed the (local) "stiffness ratio” of 
the problem. Problems may be considered to be marginally stiff if 
S(x) is 0 (10), while stiffness ratio upto 0 (10^) are not 

uncommon in practical problems arising in such fields as chemical 
kinetics, process control and reactor kinetics. 

A stiff problem is often referred to in the literature as a 
problem wilth "widely differing time constants" or as a system 
with a large "Llpscitz constant". 



1.2 GEAR’S METHOD OF SOLVING STIFF SYSTEMS 


The NAG subroutine D02EBF used for calculation Is based on 
Gear’s method. The subroutine D02EBF Is suited both for non-stiff 
as well as stiff systems- * f the required system of di f f erent ial 
equations Is non-stiff, Adams-Bashf orth p-th order predictor 
corrector technique is used. The order p of predictor - corrector 
is decided by the subroutine Itself. 


The system of N-diff erent ial equations are: 

I = My, ‘> CA.l) 

In case of Adam-Bashforth method predictor equation of order p I s 

given by 


y„ „ = y . + ft h y . + ... + ^ y n , 

n,o n-1 * 1 n 1 P n~l 


(A. 2) 


where 


y i = y< ' t I' ) = lh ^ belft 6 the ate P aiz6 ), 

» 

y L = t t L ), 

'i 


and where the ft are given, for example, in Henrici [35] - The 


approximation y in this predictor corrector scheme as the 

“» (o) 

first approximation In the Adams-Moulton corrector formula of the 

p-th order given by 


n, (m+1 ) 


y n _! o h f (y n,CnO- V * <*1 h y n-l 


+ ft . h y . .. 
p-1 n-p+1 


(A. 3) 


The coefficients ft* can be found in Henrici [35]. If the corrector 
equation (A. 3) is iterated until it converges to y ft (as is 
guaranteed for small enough atep size and smooth function O the 



1.2 GEAR* S METHOD OF SOLVING STIFF SYSTEMS 


The NAG subroutine D02EBF used for calculation is based on 
Gear’s method. The subroutine D02EBF is suited both for non-stiff 
as well as stiff systems. If the required system of differential 
equations is non-stiff, Adams-Bashf orth p-th order predictor - 
corrector technique is used. The order p of predictor - corrector 
is decided by the subroutine Itself. 


The system of N-di f f erential equations are: 
dy 


dt 


= t (y, t) 


(A . 1 ) 


In case of Adam-Bashf orth method predictor equation of order p is 
given by 


n, o 


- y „-l * ft l h y n-l + 


. + ft h y 

p n-p 


(A. 2 ) 


where 


y A = y(t £ ) = ih (h being the step size), 

» 

y i = ( (y^ t i ), 

and where the ft^ are given, for example, in Henrici [35]. The 
approximation y^ in this predictor corrector scheme as the 

first approximation in the Adams-Moulton corrector formula of the 
p-th order given by 


y , . = y „ + /3 hf(y , v , t ) + ft. hy . 
°n , (m+1 ) *n-l 'o w n,(m)' n' '1 *n-l 


+ ft . h y ... 
P-1 n-p+1 


(A. 3) 


it 

The coefficients ft. can be found in Henrici [35]. If the corrector 

l 

equation (A. 3) is iterated until it converges to y ft (as is 
guaranteed for small enough step size and smooth function f) the 



truncation error introduced in the n-th step of the integration 

will be c£ +1 h p+1 y ( ‘ p + 1) (t ft ) + 0 (h p+2 ), where y (k;) is the k-th 

A 

derivative of y. The coefficients are given in Henrlci [35]. 

The method for stiff equations is similar. It uses a p-th 
order predictor formula of the form 


n , (o ) 


ol y . + 
1 n~ 1 


+ « + n 1 h y 


n-p 


n-1 


(A. 4) 


and a corrector 


y n.c*+l) = “l y n-l + + “p y n-p + ”o h fCy n.C»0' 

(A. 5) 


The integration truncation error when (A. 4) is iterated to 

convergence is C^ +1 h P 1 y^ P 1 ^ (t^) + 0 (h P ), where C p + ^ = pTl* 
* 

The a. and r) 
l o 


are given in [38]. 



APPENDIX B 


HOPF BIFURCATION THEOREM 

Let ua consider a system o£ ordinary differential equations 
^ = f (X, «) (B.l) 

where. 


X = 

( X ^ j, 

X 2 ’ 

.... x n ) T , 


ii 


V 

.... f n ) T . 


« = 

(«j_» 

a 2 ’ 

.... « ) . 
in 

(B . 2 ) 


Here x. denote state variables, t is time, and a. is generally a 
1 

vector of system parameters. The system is autonomous, i.e., the 
time variable does not occur in the right-hand sides explicitly. 
It shall be assumed that RHS are continuous and continuously 
differentiable functions of their arguments. 

The steady-state (stationary) solution x of system (B.l) 
satisfies the set of nonlinear equations 

£ i <*!• *2 x n' = °* 1 * 1>2 " (B - 3) 

Let us consider a steady-state solution x of Eqs . (B.3) for 

« = « + and a branch of stationary solutions x(ot) in the 

neighborhood of a = « + . The functions f . are taken to be 

sufficiently smooth. 



Let ua assume that all eigen values of the Jacobian matrix J 
are nonzero and that only two eigenvalues are purely imaginary; 
let us call them X(«) and X (ot ) in the neighborhood of a = o* + . 

Then, 

Re {XO + )> = Re {X(ot + )> = 0 (B.4) 

Let us further assume that the real part of X ’ (<m + ) is non 
zero; that is. 

Re {X " (ot + ) } * 0 (B . 5 ) 

Ue then conclude that a branch of periodic solutions of (B.l) 

+ + -f 

exists for «>«,«<«, or «=**. In greater detail: then 
there exist functions X(t, *>), T(^), and et<>) defined for all t 
and for a small £, such that: 



is unstable. 
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• location of eigen values at bifurc- 
aation( arrows show movement of 
eigenvalues as a.^ increases) 


FIG 2 HOPF BIFURCATION IN 
TWO - TEMP. MODEL 
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FIG. 3 POSSIBLE ZONES OF PERIODIC 

BEHAVIOR IN TWO -TEMPERATURE 
REACTOR-MODEL WITHOUT DELAYED NEUTRONS 



Parameter values 



FIG. 4 PROJECTION OF A TRANSIENT ORBIT ON 

1T f Vs. IP PLANE WITHOUT DELAYED NEUTRANS 
(Negative fuel and coolant feedback) 



Parameter 



FIG. 5 PROJECTION OF ANOTHER TRANSIENT ORBIT ON 
TT f Vs. IP PLANE WITHOUT DELAYED NEUTRONS 
(Different values of a r and a<) 


pi a c = -50.0 p =2 

a f =-275.0 c =1/2 
9 = 1/20 



FIG. 6 PROJECTION OF A TRANSIENT ORBIT ON 
TT f Vs. IP PLANE 

(Different values of a c and a« ) 



Im (X ) 



FIG. 7 MOVEMENTS OF EIGENVALUES AS d f INCREASES 
BEYOND a* IN TWO TEMPERATURE FEEDBACK 
MODEL WITH DELAYED NEUTRONS • 

(The real eigenvalue close to the origin remains 
almost stationary in numerical calculation) 



FIG. 8 POSSIBLE ZONES OF PERIODIC BEHAVIOR IN TWO 
TEMPERATURE REACTOR MODEL WITH DELAYED 
NEUTRONS . 



= 1/20 



PROJECTION OF TRANSIENT ORBIT ON 
P Vs. IT, PLANE WITH DELAYED NEUTRONS 

(Negative fuel and coolant reactivity feedback . 
a f = -350 , a c = - 8 500). ’ 




= 1/2 
= 1/20 



FIG. 11 PROJECTION OF A TRANSIENT ORBIT 

ON TT f Vs. IP PLANE WITH DELAYED NEUTRONS 
(Negative fuel and coolant reactivity feedback $ 
a< =-240 , a r = -8000) 



= 1/20 



FIG. 12 PROJECTION OF A TRANS1ENJ pfiWJ 

ON TT f Vs. IP PLANE WITH DELAYED NEUTRONS 
(Negative fuel and coolant reactivity feedback • 
a f = -680 , a c = -10000) ' 



able 1: Summary of Time Periods and Percent Fluctuations in Power 
and Fuel Temperature for the Limit Cycles Analyzed. 

CASE 1. Without Delayed Neutrons 


xp . &c 
o . 


& £ 


Time 

Periods 

(a) 


Percent 

Power 

Fluctuation 


Percent Fuel 

Temperature 

Fluctuation 


Limit Cycle 
Plotted in 
Fig. # 


11 

6.26 

67 

o 

CO 

17 

7.46 

64 

9.0 


5 


3 . 60 


81 


5.4 


6 


3 le 2: Summary of Time Periods and Percent Fluctuations in Power 
and Fuel Temperature for the Limit Cycles Analyzed. 

CASE 2. Uith Delayed Neutrons 


a c 

a f 

Time 

Periods 

(a) 

Percent 

Power 

Fluctuation 

Percent Fuel 

Temperature 

Fluctuation 

Limit Cycle 
Plotted in 
Fifi. # 

-8500 

-350 

2.09 

400 

10.0 

9 

-9000 

-460 

1.99 

350 

9.5 

10 

-8000 

-240 

2 .23 

430 

12.0 

11 

-10000 

-680 

1.80 

300 

7.5 

12 

-6500 

89 

2.77 

570 

18.0 

13 


Table 3: Typical Values 
Present Study 
Reactors ) . 

for the Reactor Parameters Used 
(Generally Appropriate for 

In the 
Thermal 

Parameters 

Typical Values 


ft 

0.0066 


X 

0.075 s -1 


T ° 

£ 

1250 K 


^ o 

T 

c 

560 K 


T i° 

525 K 


i 

2xl0 _3 - 5xl0~ 5 a 


P o 

1000 MU 


Kl 

io" 4 - io" 6 


I a c l 

io" 3 - io‘ 6 


c. 

10 7 J k" 1 


C 

c 

2 x 10 7 J K -1 
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